using ProgressLogging using NBodySimulator, OrdinaryDiffEq, StaticArrays using Plots, DataFrames, StatsPlots function setup(t) T = 120.0 # K kb = 1.38e-23 # J/K ϵ = T * kb # J σ = 3.4e-10 # m ρ = 1374 # kg/m^3 m = 39.95 * 1.6747 * 1e-27 # kg N = 216 L = (m*N/ρ)^(1/3) R = 2.25σ v_dev = sqrt(kb * T / m) # m/s _L = L / σ _σ = 1.0 _ϵ = 1.0 _m = 1.0 _v = v_dev / sqrt(ϵ / m) _R = R / σ bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L) lj_parameters = LennardJonesParameters(_ϵ, _σ, _R) pbc = CubicPeriodicBoundaryConditions(_L) lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters)); simulation = NBodySimulation(lj_system, (0.0, t), pbc, _ϵ/T) return simulation end
setup (generic function with 1 method)
We will consider the following integrators
config(τ, at, rt) = [ # symplectic (alg=VelocityVerlet, dt=τ), (alg=McAte2, dt=τ), # (alg=CalvoSanz4, dt=τ), # (alg=McAte5, dt=τ), # (alg=Yoshida6, dt=τ), # (alg=KahanLi8, dt=τ), # DPRKN (alg=DPRKN6, abstol=at, rtol=rt), # (alg=DPRKN8, abstol=at, rtol=rt), (alg=DPRKN12, abstol=at, rtol=rt), # others (alg=Tsit5, abstol=at, rtol=rt), (alg=Vern7, abstol=at, rtol=rt), (alg=Vern9, abstol=at, rtol=rt) ]
config (generic function with 1 method)
We will first begin with the evolution of the energy error over time for different integrators. In this case we will fix both simulation time and timestep.
function benchmark(energyerr, rts, ts, t, configs) simulation = setup(t) prob = SecondOrderODEProblem(simulation) for config in configs alg = config.alg sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...) result = NBodySimulator.SimulationResult(sol, simulation) ΔE(t) = total_energy(result, t) - total_energy(result, 0) energyerr[alg] = [ΔE(t) for t in sol.t[2:end]] rts[alg] = rt ts[alg] = sol.t[2:end] end end ΔE = Dict() rt = Dict() ts = Dict() configs = config(2.3e-4, 1e-20, 1e-20) benchmark(ΔE, rt, ts, 40., configs) plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft); for c in configs plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", xscale=:log10, yscale=:log10) end plt
Now, let us consider a fixed simulation time
t = 35.0 results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[], :abstol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[]); function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs) simulation = setup(t) prob = SecondOrderODEProblem(simulation) for config in configs alg = config.alg sol, rt, b, gc, memalloc = @timed solve(prob, alg(); save_everystep=false, progress=true, progress_name="$alg", config...) result = NBodySimulator.SimulationResult(sol, simulation) ΔE = total_energy(result, t) - total_energy(result, 0) energyerr[alg] = ΔE rts[alg] = rt bytes[alg] = b allocs[alg] = memalloc nt[alg] = sol.destats.naccept nf[alg] = sol.destats.nf + sol.destats.nf2 end end
benchmark (generic function with 2 methods)
and change the timesteps. For adaptive integrators, we will use tolerances instead.
τs = 10 .^range(-4, -3, length=5) ats = 10 .^range(-20, -14, length=5) rts = 10 .^range(-20, -14, length=5)
5-element Array{Float64,1}:
1.0e-20
3.162277660168379e-19
1.0e-17
3.1622776601683793e-16
1.0e-14
The tolerances for adaptive algoritms were chosen such that the energy error is close to the one given by the VelocityVerlet with the chosen τ.
Now, let us run the benchmarks
@progress "Variable dt" for (τ, at, rt) in zip(τs, ats, rts) runtime = Dict() ΔE = Dict() nt = Dict() nf = Dict() b = Dict() allocs = Dict() GC.gc() benchmark(ΔE, runtime, b, allocs, nt, nf, t, config(τ, at, rt)) for (k,v) in ΔE push!(results, [string(k), runtime[k], at, rt, abs(v), nt[k], nf[k]]) end end
Let us now examine the results
results
| integrator | runtime | τ | abstol | EnergyError | timesteps | f_evals | |
|---|---|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Float64 | Int64 | Int64 | |
| 1 | DPRKN12 | 140.487 | 1.0e-20 | 1.0e-20 | 0.00159857 | 9729 | 177844 |
| 2 | Vern7 | 35.046 | 1.0e-20 | 1.0e-20 | 0.0688734 | 3998 | 41672 |
| 3 | DPRKN6 | 32.8794 | 1.0e-20 | 1.0e-20 | 0.0363926 | 6035 | 45238 |
| 4 | Tsit5 | 18.2855 | 1.0e-20 | 1.0e-20 | 0.1444 | 3281 | 21717 |
| 5 | VelocityVerlet | 583.327 | 1.0e-20 | 1.0e-20 | 0.00332889 | 350000 | 700002 |
| 6 | McAte2 | 588.858 | 1.0e-20 | 1.0e-20 | 0.000154077 | 350000 | 1050002 |
| 7 | Vern9 | 161.857 | 1.0e-20 | 1.0e-20 | 0.0995092 | 11191 | 191682 |
| 8 | DPRKN12 | 141.322 | 3.16228e-19 | 3.16228e-19 | 0.00443486 | 9853 | 180346 |
| 9 | Vern7 | 35.0334 | 3.16228e-19 | 3.16228e-19 | 0.111018 | 3984 | 41642 |
| 10 | DPRKN6 | 31.6658 | 3.16228e-19 | 3.16228e-19 | 0.0225537 | 5940 | 44195 |
| 11 | Tsit5 | 18.4532 | 3.16228e-19 | 3.16228e-19 | 0.141024 | 3264 | 21801 |
| 12 | VelocityVerlet | 326.416 | 3.16228e-19 | 3.16228e-19 | 0.00405704 | 196820 | 393642 |
| 13 | McAte2 | 331.795 | 3.16228e-19 | 3.16228e-19 | 0.00058266 | 196820 | 590462 |
| 14 | Vern9 | 162.485 | 3.16228e-19 | 3.16228e-19 | 0.104935 | 11173 | 191618 |
| 15 | DPRKN12 | 140.164 | 1.0e-17 | 1.0e-17 | 0.0148921 | 9737 | 177574 |
| 16 | Vern7 | 34.5065 | 1.0e-17 | 1.0e-17 | 0.0841265 | 3925 | 40712 |
| 17 | DPRKN6 | 33.174 | 1.0e-17 | 1.0e-17 | 0.0455213 | 6069 | 46043 |
| 18 | Tsit5 | 18.4812 | 1.0e-17 | 1.0e-17 | 0.304543 | 3274 | 21813 |
| 19 | VelocityVerlet | 183.214 | 1.0e-17 | 1.0e-17 | 0.00762108 | 110680 | 221362 |
| 20 | McAte2 | 184.581 | 1.0e-17 | 1.0e-17 | 0.00368266 | 110680 | 332042 |
| 21 | Vern9 | 162.971 | 1.0e-17 | 1.0e-17 | 0.0881161 | 11225 | 192226 |
| 22 | DPRKN12 | 138.089 | 3.16228e-16 | 3.16228e-16 | 0.00795424 | 9719 | 177430 |
| 23 | Vern7 | 34.945 | 3.16228e-16 | 3.16228e-16 | 0.112498 | 3959 | 41112 |
| 24 | DPRKN6 | 31.8677 | 3.16228e-16 | 3.16228e-16 | 0.0207562 | 5897 | 43922 |
| 25 | Tsit5 | 18.0651 | 3.16228e-16 | 3.16228e-16 | 0.108319 | 3249 | 21465 |
| 26 | VelocityVerlet | 104.651 | 3.16228e-16 | 3.16228e-16 | 0.00370175 | 62240 | 124482 |
| 27 | McAte2 | 103.137 | 3.16228e-16 | 3.16228e-16 | 0.00384686 | 62240 | 186722 |
| 28 | Vern9 | 160.911 | 3.16228e-16 | 3.16228e-16 | 0.0971133 | 10997 | 188226 |
| 29 | DPRKN12 | 142.737 | 1.0e-14 | 1.0e-14 | 0.0156339 | 9849 | 179842 |
| 30 | Vern7 | 34.7386 | 1.0e-14 | 1.0e-14 | 0.034158 | 3884 | 40202 |
| 31 | DPRKN6 | 31.4633 | 1.0e-14 | 1.0e-14 | 0.100083 | 5872 | 43768 |
| 32 | Tsit5 | 18.4425 | 1.0e-14 | 1.0e-14 | 0.314547 | 3259 | 21783 |
| 33 | VelocityVerlet | 58.4034 | 1.0e-14 | 1.0e-14 | 0.00366562 | 35000 | 70002 |
| 34 | McAte2 | 57.548 | 1.0e-14 | 1.0e-14 | 0.00980374 | 35000 | 105002 |
| 35 | Vern9 | 162.873 | 1.0e-14 | 1.0e-14 | 0.129173 | 11065 | 190002 |
The energy error as a function of runtime is given by
@df results plot(:EnergyError, :runtime, group=:integrator, xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")
If we consider the number of timesteps instead, we obtain
@df results plot(:EnergyError, :timesteps, group=:integrator, xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of timesteps")
and we can also see how the muber of acceleration function calls affects the runtime
@df results plot(:f_evals, :runtime, group=:integrator, xscale=:log10, yscale=:log10, xlabel="Number of f calls", ylabel="Runtime (s)")